Effects of Bond Yields, VIX and Initial Jobless Claims on S&P 500

Updated on: 2022-07-30

1 Introduction

The purpose of this project is to study the combined effects of bond yields, in particular the 10Y/2Y Treasury Yield Spread, the CBOE Volatility Index, and the U.S. Initial Jobless Claims on the S&P 500, measured using the SPDR S&P 500 ETF (SPY) closing price. It uses an Autoregressive Distributed Lags model to estimate the SPY returns and forecast the returns for the 4 weeks in June 2022.

2 Packages Required

library(ARDL) # For ARDL models lag selection
library(corrplot) # For visualizing correlation between variables
library(dLagM) # For ARDL models and forecasting
library(dplyr) # For data manipulation
library(forecast) # For ARIMA model and evaluation metrics
library(lubridate) # For manipulating time series objects
library(quantmod) # For obtaining historical data from Yahoo Finance and FRED
library(urca) # For unit root tests
library(stargazer) # For tidy regression output and tables, where possible

3 Description of Variables

In this section, I briefly described the four variables that were used in the project and their expected univariate relationship with the S&P 500. I would be using the weekly data of each variable since data for Initial Jobless Claims is released weekly.

3.1 SPDR S&P 500 ETF

The S&P 500 Index (ticker in Yahoo Finance: ^GSPC) measures the stock performance of the 500 largest companies listed in the U.S. by market capitalization and is usually used as a gauge of the overall U.S. stock market. The SPDR S&P 500 ETF (ticker: SPY) tracks this index and will be used as the dependent variable in the models that will be used in Section 4.

I have downloaded and stored the data into the object spy_price.

# Retrieve SPY historical data from Yahoo Finance using quantmod package

# Set start and end date for data retrieval
startdate <- as.Date("2000-01-03") # First trading day of year 2000
enddate <- as.Date("2022-07-01") 

spy_price <- getSymbols(Symbols = "SPY", 
                        src = "yahoo", 
                        auto.assign = F, 
                        from = startdate, to = enddate, 
                        periodicity = "weekly")

# Show first and last 6 observations
head(spy_price); tail(spy_price)

##            SPY.Open SPY.High  SPY.Low SPY.Close SPY.Volume SPY.Adjusted
## 2000-01-03 148.2500 148.2500 137.2500  145.7500   42725700     96.34634
## 2000-01-10 146.2500 147.4688 142.8750  146.9688   32748700     97.15200
## 2000-01-17 145.3438 147.0000 143.8125  144.4375   24691300     95.47875
## 2000-01-24 145.6562 145.8438 135.5312  135.8750   45836400     89.81861
## 2000-01-31 135.8125 144.0000 135.0000  142.5938   38317400     94.25992
## 2000-02-07 142.5625 144.5625 138.0312  138.6875   35834100     91.67777

##            SPY.Open SPY.High SPY.Low SPY.Close SPY.Volume SPY.Adjusted
## 2022-05-23   392.83   415.38  386.96    415.26  426273600     413.4739
## 2022-05-30   413.55   417.44  406.93    410.54  334006700     408.7742
## 2022-06-06   414.78   416.61  389.75    389.80  400315000     388.1234
## 2022-06-13   379.85   383.90  362.17    365.86  645270700     364.2864
## 2022-06-20   371.89   390.09  370.18    390.08  344213700     390.0800
## 2022-06-27   391.05   393.16  372.56    377.25  330742800     377.2500

# Number of rows of data
nrow(spy_price)

## [1] 1174

# Check for missing data
colSums(is.na(spy_price))

##     SPY.Open     SPY.High      SPY.Low    SPY.Close   SPY.Volume SPY.Adjusted 
##            0            0            0            0            0            0

The dates can be a little misleading when downloading weekly data but here is what each column refers to:

Date: Date of the observation, set as Monday (first trading day of year 2000)
SPY.Open: Opening price of the week on Monday
SPY.High: Highest price during the week
SPY.Low: Lowest price during the week
SPY.Close: Closing price of the week on Friday (except last observation which is last trading day of June 2022), adjusted for splits
SPY.Volume: Total volume traded during the week
SPY.Adjusted: Closing price of the week (except last observation which is last trading day of June 2022), adjusted for splits and dividends

If we want returns to be similar to the S&P 500 Index, the unadjusted closing price should be used since the index does not give the total return which includes dividends. Accounting for dividends would reduce losses during bad times and increase gains during good times, which does not reflect the week-to-week movement of the index. Since the purpose of the project is to predict U.S. stock market, unadjusted closing price is used.

Plot of SPY over the years:

spy_close <- spy_price[,"SPY.Close"] %>% 
  `colnames<-`("SPY")

plot(spy_close, main = "SPY Weekly Closing Price")

Since this is a time series analysis, we should also check for stationarity of our data. The chart shows that the series is not stationary and it can be checked with the Augmented Dickey Fuller and KPSS Tests.

# Add deterministic trend to tests since the chart seems to show it

spy_close %>% urca::ur.df(type = "trend", selectlags = "AIC") %>% summary()

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -43.320  -1.818   0.212   2.086  28.827 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  0.0348264  0.3180695   0.109  0.91283   
## z.lag.1     -0.0053702  0.0028262  -1.900  0.05766 . 
## tt           0.0019720  0.0007691   2.564  0.01047 * 
## z.diff.lag  -0.0758051  0.0292449  -2.592  0.00966 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.89 on 1168 degrees of freedom
## Multiple R-squared:  0.01139,    Adjusted R-squared:  0.008849 
## F-statistic: 4.485 on 3 and 1168 DF,  p-value: 0.003877
## 
## 
## Value of test-statistic is: -1.9002 2.9906 3.3841 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.96 -3.41 -3.12
## phi2  6.09  4.68  4.03
## phi3  8.27  6.25  5.34

spy_close %>% urca::ur.kpss(type = "tau", lags = "long") %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 22 lags. 
## 
## Value of test-statistic is: 1.0198 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216

ADF test is unable to reject the null hypothesis which claimed that there is unit root (-1.90 vs -3.41, left-tailed test) and KPSS test rejected the null hypothesis which claimed that it is stationary (1.02 vs 0.146). I did the tests again with the weekly returns calculated from the unadjusted closing price using the simple/discrete method.

# Calculate discrete returns, na.omit to remove first observation as it will return NA
returns_spy <- na.omit(diff(x = spy_close, lag = 1, differences = 1) / stats::lag(x = spy_close, k = 1))

nrow(returns_spy)

## [1] 1173

plot(returns_spy, main = "SPY Weekly Returns")

returns_spy %>% urca::ur.df(type = "drift", selectlags = "AIC") %>% summary()

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.205086 -0.011736  0.001236  0.013325  0.127818 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0011954  0.0007351   1.626    0.104    
## z.lag.1     -1.0457649  0.0430236 -24.307   <2e-16 ***
## z.diff.lag  -0.0326059  0.0293380  -1.111    0.267    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0251 on 1168 degrees of freedom
## Multiple R-squared:  0.5407, Adjusted R-squared:  0.5399 
## F-statistic: 687.4 on 2 and 1168 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -24.3068 295.4093 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78

returns_spy %>% urca::ur.kpss(type = "mu", lags = "long") %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 22 lags. 
## 
## Value of test-statistic is: 0.3113 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

With the discrete returns, ADF test rejects the null hypothesis (-24.31 vs -2.86) and KPSS test is unable to reject the null (0.31 vs 0.463). The return series, which is simply the first difference of the price data, is stationary.

3.2 10Y/2Y Treasury Yield Spread

The yield spread is the difference between the 2-year and 10-year treasury bonds, and is watched by many professionals as an early indicator of stock market downturns when the spread turns negative. This is because longer-dated bonds should have higher yields than bonds with shorter durations, so it is normal when the 10-year treasury yield is above the 2-year treasury yield. However, if the 2-year yield is higher than the 10-year yield, we have an inverted yield curve which signals that investors expect long-term interest rates to fall. This means that the yield spread and the S&P 500 should be negatively correlated.

I have downloaded and stored the yield spread data into the object tys.

# Retrieve 10Y/2Y Treasury Yield Spread from St. Louis Fed's FRED using quantmod package

tys <- getSymbols(Symbols = "T10Y2Y", src = "FRED", auto.assign = F, from = startdate, to = enddate, periodicity = "weekly")

head(tys); tail(tys)

##            T10Y2Y
## 1976-06-01   0.68
## 1976-06-02   0.71
## 1976-06-03   0.70
## 1976-06-04   0.77
## 1976-06-07   0.79
## 1976-06-08   0.79

##            T10Y2Y
## 2022-07-22  -0.21
## 2022-07-25  -0.19
## 2022-07-26  -0.21
## 2022-07-27  -0.18
## 2022-07-28  -0.17
## 2022-07-29  -0.22

The data retrieved from FRED was not automatically converted by the quantmod package to weekly data, as compared to the SPY ETF data retrieved from Yahoo Finance. Furthermore, the from and to arguments does not work when retrieving FRED data. I manually adjusted this to match the SPY price data output above.

# 1. Adjust data collected to start from 2000 and end in 2022

tys <- tys["2000/2022-06",] %>% 
  `colnames<-`("TYS")

# 2. Check for missing data. Replace NAs with prior observation.

sum(is.na(tys))

## [1] 241

clean_tys <- na.locf(object = tys)

# 3. Change daily frequency to weekly frequency (since SPY closing price is Friday, we extract Friday data)
# But last observation in closing price is the last trading day of June 2022, so remember to add that in

# Indicate week_start = 1 for week to start on Monday
weekly_tys <- rbind(clean_tys[wday(clean_tys, week_start = 1) == 5], last(clean_tys))

# Third, adjust dates since second step will return dates on Friday

index(weekly_tys) <- index(spy_close)

head(weekly_tys); tail(weekly_tys); nrow(weekly_tys)

##              TYS
## 2000-01-03  0.21
## 2000-01-10  0.25
## 2000-01-17  0.31
## 2000-01-24  0.08
## 2000-01-31 -0.10
## 2000-02-07 -0.02

##             TYS
## 2022-05-23 0.27
## 2022-05-30 0.30
## 2022-06-06 0.09
## 2022-06-13 0.08
## 2022-06-20 0.09
## 2022-06-27 0.06

## [1] 1174

Plot SPY Closing Price and 10Y/2Y Treasury Yield Spread over time:

plot(merge(spy_close, weekly_tys), multi.panel = T, yaxis.same = F, main = "SPY Closing Price and 10Y/2Y Treasury Yield Spread")

We can see that the yield spread is negative before the stock market collapsed in the early 2000s and in 2008/2009, and the spread peaked around the bottom of the stock market. The yield spread does not seem to be stationary.

weekly_tys %>% urca::ur.df(type = "drift", selectlags = "AIC") %>% summary()

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.42238 -0.04675 -0.00585  0.04415  0.43415 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  0.004788   0.004245   1.128   0.2596  
## z.lag.1     -0.003959   0.002764  -1.432   0.1524  
## z.diff.lag  -0.065091   0.029190  -2.230   0.0259 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08413 on 1169 degrees of freedom
## Multiple R-squared:  0.00625,    Adjusted R-squared:  0.00455 
## F-statistic: 3.676 on 2 and 1169 DF,  p-value: 0.02561
## 
## 
## Value of test-statistic is: -1.432 1.0276 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78

weekly_tys %>% urca::ur.kpss(type = "mu", lags = "long") %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 22 lags. 
## 
## Value of test-statistic is: 0.5241 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

ADF test shows that we cannot reject the null hypothesis (-1.43 vs -2.86), so unit root is present in the series. KPSS test shows that the null hypothesis is rejected (0.52 vs 0.463), so the series is not stationary.

Since the yield spread is stated in percentage, I would simply take the first difference which would refer to the week-to-week changes in yield spread.

# Remove first observation since it will return NA after taking first difference
diff_tys <- na.omit(diff(x = weekly_tys, lag = 1, differences = 1))

nrow(diff_tys)

## [1] 1173

plot(diff_tys, main = "First-Difference of 10Y/2Y Treasury Yield Spread")

diff_tys %>% urca::ur.df(type = "drift", selectlags = "AIC") %>% summary()

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41751 -0.04634 -0.00487  0.04433  0.43262 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0002182  0.0024595  -0.089    0.929    
## z.lag.1     -1.0276976  0.0427065 -24.064   <2e-16 ***
## z.diff.lag  -0.0372397  0.0292319  -1.274    0.203    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08416 on 1168 degrees of freedom
## Multiple R-squared:  0.5344, Adjusted R-squared:  0.5336 
## F-statistic: 670.4 on 2 and 1168 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -24.0642 289.5428 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78

diff_tys %>% urca::ur.kpss(type = "mu", lags = "long") %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 22 lags. 
## 
## Value of test-statistic is: 0.185 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

After taking first difference, the number of observations is the same as returns_spy, and ADF and KPSS test show that the series is stationary.

3.3 CBOE Volatility Index

The CBOE Volatility Index, commonly known as VIX, is used to assess volatility of the S&P 500 Index. A higher VIX value would mean greater fear and uncertainty in the market.

The VIX data is stored into vix.

# Retrieve VIX from FRED, adjustment of data is similar to the yield spread

vix <- getSymbols(Symbols = "VIXCLS", src = "FRED", auto.assign = F)

vix <- vix["2000/2022-06",] %>%
  `colnames<-`("VIX")

# Check for missing data and replace with prior observation
sum(is.na(vix))

## [1] 206

clean_vix <- na.locf(object = vix)

weekly_vix <- rbind(clean_vix[wday(clean_vix, week_start = 1) == 5], last(clean_vix))

index(weekly_vix) <- index(spy_close)

head(weekly_vix); tail(weekly_vix); nrow(weekly_vix)

##              VIX
## 2000-01-03 21.72
## 2000-01-10 19.66
## 2000-01-17 20.82
## 2000-01-24 26.14
## 2000-01-31 21.54
## 2000-02-07 24.42

##              VIX
## 2022-05-23 25.72
## 2022-05-30 24.79
## 2022-06-06 27.75
## 2022-06-13 31.13
## 2022-06-20 27.23
## 2022-06-27 28.71

## [1] 1174

Plot SPY Closing Price and VIX over time:

plot(merge(spy_close, weekly_vix), multi.panel = T, yaxis.same = F, main = "SPY Closing Price and VIX")

The chart shows that VIX is higher during stock market crashes, and tends to be around 10 to 30 when the stock market is growing (2004-2007 and 2013-2019). The correlation between VIX and S&P 500 is expected to be negative.

Check for stationarity of VIX series:

weekly_vix %>% urca::ur.df(type = "drift", selectlags = "AIC") %>% summary()

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.8603  -1.5482  -0.3863   1.1857  27.5899 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.28935    0.24050   5.361 9.95e-08 ***
## z.lag.1     -0.06446    0.01112  -5.799 8.58e-09 ***
## z.diff.lag  -0.11149    0.02907  -3.835 0.000133 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.256 on 1169 degrees of freedom
## Multiple R-squared:  0.04808,    Adjusted R-squared:  0.04645 
## F-statistic: 29.52 on 2 and 1169 DF,  p-value: 3.11e-13
## 
## 
## Value of test-statistic is: -5.7989 16.8176 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78

weekly_vix %>% urca::ur.kpss(type = "mu", lags = "long") %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 22 lags. 
## 
## Value of test-statistic is: 0.3443 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

Based on the ADF and KPSS tests, the VIX series at level is stationary with a drift.

3.4 U.S. Initial Jobless Claims

The U.S. Initial Jobless Claims is a weekly economic data measuring the number of people filing for unemployment claims for the first time in the past week. A growth in initial jobless claims is typically seen just before and during a recession, and gradually decline as the economy recovers.

I stored the data into ijc.

# Retrieve initial jobless claims from FRED, adjustment of data is similar to the yield spread

ijc <- getSymbols(Symbols = "ICSA", src = "FRED", auto.assign = F)

ijc <- ijc["2000/2022-06",] %>%
  `colnames<-`("IJC")

# Check for missing data
sum(is.na(ijc))

## [1] 0

nrow(ijc)

## [1] 1174

# Initial Jobless Claims is weekly data, but adjust the dates so that it can be plotted
# and merged into same object with other data later.
index(ijc) <- index(spy_close)

head(ijc); tail(ijc)

##               IJC
## 2000-01-03 286000
## 2000-01-10 298000
## 2000-01-17 289000
## 2000-01-24 284000
## 2000-01-31 285000
## 2000-02-07 312000

##               IJC
## 2022-05-23 211000
## 2022-05-30 202000
## 2022-06-06 232000
## 2022-06-13 231000
## 2022-06-20 233000
## 2022-06-27 231000

Plot SPY Closing Price and Initial Jobless Claims over time:

plot(merge(spy_close, ijc), multi.panel = T, yaxis.same = F, main = "SPY Closing Price and Initial Jobless Claims")

Due to the COVID-19 pandemic beginning in early 2020, there was an influx of initial unemployment claims that caused the historical patterns of claims to look flat. Such an increase in numbers may affect the estimated models.

plot(merge(spy_close, ijc)["/2019",], multi.panel = T, yaxis.same = F, main = "SPY and Initial Jobless Claims Before 2020")

The chart of SPY and Initial Jobless Claims before 2020 shows unemployment claims is high during the stock market crashes in early 2000s and in 2008/2009. This makes sense, and we expect the unemployment claims to be negatively correlated with the S&P 500.

To check for stationarity in the series, I only used data before 2020 since the extremely high unemployment claims during the COVID-19 pandemic will likely cause the series to be tested as stationary even when it is not.

ijc["/2019",] %>% urca::ur.df(type = "drift", selectlags = "AIC") %>% summary()

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -74891  -7982   -472   7554  97164 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.874e+03  1.863e+03   2.079   0.0378 *  
## z.lag.1     -1.138e-02  5.211e-03  -2.185   0.0291 *  
## z.diff.lag  -1.658e-01  3.060e-02  -5.417 7.51e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14940 on 1039 degrees of freedom
## Multiple R-squared:  0.03394,    Adjusted R-squared:  0.03208 
## F-statistic: 18.25 on 2 and 1039 DF,  p-value: 1.623e-08
## 
## 
## Value of test-statistic is: -2.1847 2.3977 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78

ijc["/2019",] %>% urca::ur.kpss(type = "mu", lags = "long") %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 21 lags. 
## 
## Value of test-statistic is: 1.51 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

The ADF test tells us that the null hypothesis cannot be rejected (-2.18 vs -2.86) and KPSS test tells us that the null hypothesis is rejected (1.51 vs 0.463). Therefore, the series is not stationary at levels.

I conducted the ADF and KPSS tests on the growth in initial jobless claims, which is the first-difference of the series.

# Take difference of data at time t and t-1 and divide by data at time t-1
# Gives us the growth in initial jobless claims (in decimal format)
# Remove first observation since it will return NA 
ijc_growth <- na.omit(diff(x = ijc, lag = 1, differences = 1) / stats::lag(x = ijc, k = 1))

head(ijc_growth)

##                     IJC
## 2000-01-10  0.041958042
## 2000-01-17 -0.030201342
## 2000-01-24 -0.017301038
## 2000-01-31  0.003521127
## 2000-02-07  0.094736842
## 2000-02-14 -0.038461538

nrow(ijc_growth)

## [1] 1173

ijc_growth["/2019",] %>% urca::ur.df(type = "drift", selectlags = "AIC") %>% summary()

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.131549 -0.024894 -0.002861  0.022950  0.301255 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.001031   0.001268   0.813    0.416    
## z.lag.1     -1.469388   0.045889 -32.021   <2e-16 ***
## z.diff.lag   0.255909   0.029986   8.534   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0409 on 1038 degrees of freedom
## Multiple R-squared:  0.6122, Adjusted R-squared:  0.6115 
## F-statistic: 819.5 on 2 and 1038 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -32.0205 512.6579 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78

ijc_growth["/2019",] %>% urca::ur.kpss(type = "mu", lags = "long") %>% summary()

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 21 lags. 
## 
## Value of test-statistic is: 0.2218 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

The ADF and KPSS tests show that the growth rate in initial jobless claims is stationary.

3.5 Summary

To summarize what had been done in this section, I discussed how the yield spread, the VIX and initial jobless claims would correlate with the S&P 500. I also carried out the ADF and KPSS tests and found that only the VIX is stationary at levels, while the other three variables are stationary at first difference. Understanding the order of integration allows us to find out if we should test for cointegration, which is not needed unless we are are only estimating the effects of the yield spread and initial jobless claims on the S&P 500.

The correlation heatmaps shows how the data at levels and the stationary data correlate with other variables in their respective groups.

corr_level <- cor(x = merge(spy_close, weekly_tys, weekly_vix, ijc), method = "spearman")

corrplot(corr = corr_level, method = "color", type = "lower", title = "Correlation of variables at level", addCoef.col = "black", mar = c(1,1,2,1))

# Remove first row of weekly_vix because it did not require differencing, so it had an additional first observation
corr_stationary <- cor(x = merge(returns_spy, diff_tys, weekly_vix[-1,], ijc_growth), method = "spearman")

corrplot(corr = corr_stationary, method = "color", type = "lower", title = "Correlation of stationary variables", addCoef.col = "black", mar = c(1,1,2,1))

When all variables are at levels, TYS, VIX and IJC were negatively correlated to SPY as I had expected. But the stationary variables paint a different picture. The weekly change in TYS was positively correlated with SPY returns, VIX negatively correlated with SPY returns, and weekly percentage change in IJC was slightly negatively correlated with SPY returns.

4 Building the Regression Models

In this section, three different types of regression models will be created with returns_spy as the dependent variable, namely an ARMA model, ARDL models with one independent variable and ARDL model with all independent variables.

For all model types, I used stationary data up to end-May 2022 as a training set and June 2022 as the test set to evaluate forecast performance.

4.1 Autoregressive Moving Average (ARMA) Model

The ARMA model is the simplest model to estimate, since it only uses the lags of the dependent variable and the lags of the error term. A plot of the autocorrelation function (ACF) and partial autocorrelation function (PACF) of the dependent variable can provide some information about the number of lags to choose from.

par(mfrow = c(2,1), mar = c(2,3,4,2))

forecast::Acf(x = returns_spy["/2022-05",], main = "ACF of SPY Returns")

forecast::Pacf(x = returns_spy["/2022-05",], main = "PACF of SPY Returns")

The ACF and PACF suggested an ARMA(3, 3) or ARMA(4, 4) model would be a good starting point to select the number of AR and MA lags, but that seems to be a very long model. Therefore, I will simply use the auto.arima function to choose an ARMA(p, q) model with the lowest Akaike Information Criterion (AIC).

arma <- forecast::auto.arima(y = returns_spy["/2022-05",], 
                             max.p = 4, max.q = 4, 
                             ic = "aic", 
                             stepwise = F, approximation = F, 
                             trace = F)

arma

## Series: returns_spy["/2022-05", ] 
## ARIMA(0,0,4) with non-zero mean 
## 
## Coefficients:
##           ma1     ma2      ma3      ma4    mean
##       -0.0764  0.0463  -0.0762  -0.0579  0.0012
## s.e.   0.0292  0.0294   0.0292   0.0288  0.0006
## 
## sigma^2 = 0.000617:  log likelihood = 2663.6
## AIC=-5315.19   AICc=-5315.12   BIC=-5284.81

The estimated ARMA(p, q) model is:

\[ \hat{r}_t = 0.0012 - 0.0764u_{t-1} + 0.0463u_{t-2} - 0.0762u_{t-3} - 0.0579u_{t-4} \]

Plot of actual and fitted values:

plot.zoo(returns_spy["/2022-05",], 
         col = "black", type = "l", lwd = 2, 
         ylab = "SPY Returns", xlab = "Time", 
         main = "Fitted SPY Returns Using ARMA")

lines(zoo(x = fitted(arma), order.by = index(returns_spy["/2022-05",])), col = "red", type = "l", lwd = 1.5)

legend(x = "bottomleft", legend = c("Actual Returns", "Fitted Values"), col = c("black", "red"), lwd = 2)

4.2 Autoregressive Distributed Lag (ARDL) Model With One Independent Variable

I would build three ARDL models, one for each independent variable, to test how the variables affect the SPY returns individually. For this, I would use the auto_ardl function in the ARDL package.

4.2.1 Regress SPY Returns on Treasury Yield Spread

# max_order = 4 to indicate maximum lag in the search for all variables is 4
# possible to state a vector of length equal to no. of variables in max_order
# grid = T to to prevent stepwise search of models

spytys <- ARDL::auto_ardl(formula = SPY ~ TYS, 
                          data = as.zoo(merge(returns_spy, diff_tys)["/2022-05",]), 
                          max_order = 4, 
                          selection = "AIC", grid = T)

summary(spytys$best_model)

## 
## Time series regression with "zoo" data:
## Start = 2000-01-17, End = 2022-05-30
## 
## Call:
## dynlm::dynlm(formula = full_formula, data = data, start = start, 
##     end = end)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.201475 -0.011267  0.001883  0.013216  0.118899 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0012728  0.0007247   1.756  0.07931 .  
## L(SPY, 1)   -0.0625140  0.0291856  -2.142  0.03240 *  
## TYS          0.0235748  0.0086463   2.727  0.00650 ** 
## L(TYS, 1)   -0.0304350  0.0086484  -3.519  0.00045 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02474 on 1164 degrees of freedom
## Multiple R-squared:  0.02388,    Adjusted R-squared:  0.02136 
## F-statistic: 9.491 on 3 and 1164 DF,  p-value: 3.389e-06

The estimated model is:

\[ \hat{r}_t = 0.0013 - 0.0625r_{t-1} + 0.0236TYS_t - 0.0304TYS_{t-1} \]

Based on the estimated coefficients, the Treasury Yield Spread at time \(t\) has positive impact on SPY returns but the Treasury Yield Spread at time \(t-1\) has a negative impact on SPY returns. The low Adjusted R-squared value means that the weekly change in TYS does not explain SPY returns well.

Plot of actual and fitted values:

plot.zoo(returns_spy["/2022-05",], 
         col = "black", type = "l", lwd = 2, 
         ylab = "SPY Returns", xlab = "Time", 
         main = "Regressing SPY Returns on Weekly Change in TYS")

# NA for first value as the lag used reduces observations by 1
lines(fitted(spytys$best_model), col = "red", type = "l", lwd = 1.5)

legend(x = "bottomleft", legend = c("Actual Returns", "Fitted Values"), col = c("black", "red"), lwd = 2)

4.2.2 Regress SPY Returns on VIX

spyvix <- ARDL::auto_ardl(formula = SPY ~ VIX, 
                          data = as.zoo(merge(returns_spy, weekly_vix[-1,])["/2022-05",]), 
                          max_order = 4, 
                          selection = "AIC", grid = T)

summary(spyvix$best_model)

## 
## Time series regression with "zoo" data:
## Start = 2000-02-07, End = 2022-05-30
## 
## Call:
## dynlm::dynlm(formula = full_formula, data = data, start = start, 
##     end = end)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.100654 -0.007445 -0.000257  0.006955  0.113733 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0077923  0.0012223   6.375 2.64e-10 ***
## L(SPY, 1)   -0.1383703  0.0291963  -4.739 2.41e-06 ***
## L(SPY, 2)   -0.0096003  0.0294643  -0.326 0.744614    
## L(SPY, 3)   -0.1013077  0.0291234  -3.479 0.000523 ***
## VIX         -0.0060032  0.0001371 -43.800  < 2e-16 ***
## L(VIX, 1)    0.0044331  0.0002491  17.797  < 2e-16 ***
## L(VIX, 2)    0.0009027  0.0002819   3.202 0.001403 ** 
## L(VIX, 3)   -0.0004992  0.0002816  -1.773 0.076533 .  
## L(VIX, 4)    0.0008522  0.0002162   3.941 8.59e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01518 on 1156 degrees of freedom
## Multiple R-squared:  0.6318, Adjusted R-squared:  0.6293 
## F-statistic:   248 on 8 and 1156 DF,  p-value: < 2.2e-16

The estimated model is:

\[ \begin{aligned} \hat{r}_t = &0.0078 -0.1384r_{t-1} -0.0096r_{t-2} -0.1013r_{t-3} \\ &-0.0060VIX_t + 0.0044VIX_{t-1} + 0.0009VIX_{t-2} -0.0005VIX_{t-3} + 0.0008VIX_{t-4} \end{aligned} \]

Unlike what we had thought about the negative correlation between SPY returns and VIX, most of the lags of VIX suggested that higher VIX values generally increases SPY returns.

Plot of actual and fitted values:

plot.zoo(returns_spy["/2022-05",], 
         col = "black", type = "l", lwd = 2, 
         ylab = "SPY Returns", xlab = "Time", 
         main = "Regressing SPY Returns on VIX")

# NA for first 4 values as the lag used reduces observations by 4
lines(fitted(spyvix$best_model), col = "red", type = "l", lwd = 1.5)

legend(x = "bottomleft", legend = c("Actual Returns", "Fitted Values"), col = c("black", "red"), lwd = 2)

4.2.3 Regress SPY Returns on Initial Jobless Claims

spyijc <- ARDL::auto_ardl(formula = SPY ~ IJC, 
                          data = as.zoo(merge(returns_spy, ijc_growth)["/2022-05",]), 
                          max_order = 4, 
                          selection = "AIC", grid = T)

summary(spyijc$best_model)

## 
## Time series regression with "zoo" data:
## Start = 2000-01-31, End = 2022-05-30
## 
## Call:
## dynlm::dynlm(formula = full_formula, data = data, start = start, 
##     end = end)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.204112 -0.011094  0.001525  0.012876  0.128448 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0012073  0.0007229   1.670 0.095198 .  
## L(SPY, 1)   -0.0523752  0.0295871  -1.770 0.076956 .  
## L(SPY, 2)    0.0326901  0.0298948   1.094 0.274400    
## L(SPY, 3)   -0.0410731  0.0296089  -1.387 0.165651    
## IJC          0.0075333  0.0020526   3.670 0.000253 ***
## L(IJC, 1)   -0.0020063  0.0020737  -0.967 0.333498    
## L(IJC, 2)    0.0085988  0.0020423   4.210 2.75e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.02455 on 1159 degrees of freedom
## Multiple R-squared:  0.0375, Adjusted R-squared:  0.03251 
## F-statistic: 7.525 on 6 and 1159 DF,  p-value: 6.245e-08

The estimated model is:

\[ \hat{r}_t = 0.0012 -0.0524r_{t-1} + 0.0327r_{t-2} -0.0411r_{t-3} + 0.0075IJC_t -0.0020IJC_{t-1} + 0.0086IJC_{t-2} \]

The low Adjusted R-Squared value means that the weekly percentage change in Initial Jobless Claims does not explain SPY Returns well. If we simply look at the coefficients of IJC and its lags, an increase in initial jobless claims seemed to generally increase SPY returns.

Plot of actual and fitted values:

plot.zoo(returns_spy["/2022-05",], 
         col = "black", type = "l", lwd = 2, 
         ylab = "SPY Returns", xlab = "Time", 
         main = "Regressing SPY Returns on Weekly Percentage Change in IJC")

# NA for first 3 values as the lag used reduces observations by 3
lines(fitted(spyijc$best_model), col = "red", type = "l", lwd = 1.5)

legend(x = "bottomleft", legend = c("Actual Returns", "Fitted Values"), col = c("black", "red"), lwd = 2)

4.3 ARDL Model With All Independent Variables

The process for building this model is the same as the ARDL models with one independent variable.

spyall <- ARDL::auto_ardl(formula = SPY ~ TYS + VIX + IJC, 
                          data = as.zoo(merge(returns_spy, diff_tys, weekly_vix[-1,], ijc_growth)["/2022-05",]), 
                          max_order = 4, 
                          selection = "AIC", grid = T)

summary(spyall$best_model)

## 
## Time series regression with "zoo" data:
## Start = 2000-02-07, End = 2022-05-30
## 
## Call:
## dynlm::dynlm(formula = full_formula, data = data, start = start, 
##     end = end)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.073324 -0.007520 -0.000163  0.006493  0.082511 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.665e-03  1.191e-03   6.436 1.80e-10 ***
## L(SPY, 1)   -2.092e-02  2.962e-02  -0.706 0.480217    
## L(SPY, 2)   -5.007e-02  2.962e-02  -1.691 0.091185 .  
## L(SPY, 3)   -1.114e-01  2.907e-02  -3.832 0.000134 ***
## L(SPY, 4)    2.573e-02  1.815e-02   1.418 0.156499    
## TYS          1.434e-02  5.050e-03   2.840 0.004591 ** 
## VIX         -6.152e-03  1.304e-04 -47.168  < 2e-16 ***
## L(VIX, 1)    5.217e-03  2.500e-04  20.867  < 2e-16 ***
## L(VIX, 2)    2.772e-05  2.971e-04   0.093 0.925678    
## L(VIX, 3)   -3.771e-04  2.958e-04  -1.275 0.202650    
## L(VIX, 4)    9.686e-04  2.180e-04   4.442 9.76e-06 ***
## IJC          1.032e-02  1.220e-03   8.460  < 2e-16 ***
## L(IJC, 1)   -1.061e-02  1.262e-03  -8.402  < 2e-16 ***
## L(IJC, 2)    7.210e-03  1.302e-03   5.536 3.82e-08 ***
## L(IJC, 3)    6.406e-04  1.306e-03   0.490 0.623885    
## L(IJC, 4)   -3.706e-03  1.268e-03  -2.923 0.003537 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01423 on 1149 degrees of freedom
## Multiple R-squared:  0.6785, Adjusted R-squared:  0.6743 
## F-statistic: 161.7 on 15 and 1149 DF,  p-value: < 2.2e-16

Based on ARDL automatic selection, the model with the lowest AIC is rather complex. In short, adding 4 lags each of SPY returns, VIX and initial jobless claims, as well as the contemporaneous treasury yield spread leads to the best model.

Plot of actual and fitted values:

plot.zoo(returns_spy["/2022-05",], 
         col = "black", type = "l", lwd = 2, 
         ylab = "SPY Returns", xlab = "Time", 
         main = "Regressing SPY Returns on All Independent Variables")

# NA for first 4 values as the lag used reduces observations by 4
lines(fitted(spyall$best_model), col = "red", type = "l", lwd = 1.5)

legend(x = "bottomleft", legend = c("Actual Returns", "Fitted Values"), col = c("black", "red"), lwd = 2)

5 Model Evaluation

In this section, I evaluated the out-of-sample forecasting accuracy for the ARMA model and the ARDL with all independent variables model. I forecasted 4 periods ahead since the test set containing June 2022 SPY returns only has 4 observations.

5.1 ARMA Model

f_arma <- forecast::forecast(object = arma, h = 4, level = 95)

stargazer(as.data.frame(f_arma), 
          type = "text", 
          title = "Forecasted SPY Returns for June 2022 Using ARMA Model",
          summary = F)

## 
## Forecasted SPY Returns for June 2022 Using ARMA Model
## ================================
##      Point Forecast Lo 95  Hi 95
## --------------------------------
## 8184     0.009      -0.040 0.058
## 8191     -0.002     -0.050 0.047
## 8198     -0.002     -0.051 0.047
## 8205     0.002      -0.047 0.051
## --------------------------------

stargazer(forecast::accuracy(f_arma, x = returns_spy["2022-06"]), 
          type = "text", 
          title = "Evaluation of ARMA Model Forecast")

## 
## Evaluation of ARMA Model Forecast
## ==============================================================
##                 ME    RMSE   MAE    MPE    MAPE   MASE   ACF1 
## --------------------------------------------------------------
## Training set -0.00000 0.025 0.017  -Inf     Inf   0.993 -0.001
## Test set      -0.022  0.057 0.055 105.732 105.732 3.204       
## --------------------------------------------------------------

For the evaluation metrics, Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE) is the simplest to understand. The Mean Percentage Error (MPE) and Mean Absolute Percentage Error (MAPE) is not particularly useful in our case since the forecasted values tend to be very small, which inflates these measures.

Plot forecasted and actual SPY at level (price data):

# Add in-sample periods Jan to May 2022, and out-of-sample period Jun 2022 (26 observations)
# Bind actual with forecasted values 

f_arma.level <- zoo(x = cbind(fspy = as.vector(last(spy_close["2022-05",])) * cumprod(1 + f_arma$mean),
                              fspy_lower = as.vector(last(spy_close["2022-05",])) * cumprod(1 + f_arma$lower),
                              fspy_upper = as.vector(last(spy_close["2022-05",])) * cumprod(1 + f_arma$upper)),
                    order.by = index(spy_close["2022-06",]))

plot.zoo(cbind(f_arma.level, spy_close["2022-06",], spy_close["2022-01/2022-05",]), 
         plot.type = "single", 
         col = c("red", "gray", "gray", "blue", "black"), lwd = 2, 
         ylab = "SPY Weekly Close Price", main = "Actual and Forecasted SPY Price Using ARMA Model")

legend(x = "bottomleft", 
       legend = c("Forecasted Price", "Lower and Upper Bound", "Actual Price June 2022", "Actual Price in Train Period"), 
       col = c("red", "gray", "blue", "black"), lwd = 1.5)

5.2 ARDL Model

5.2.1 Re-estimating ARDL Model

For the ARDL model, more steps are required to obtain a forecast since the forecast function used in the ARMA model does not work with ARDL models, and the ARDL package does not have a function for it. I would be using the dLagM package and the forecast function available in it.

# Using dLagM package to create same ARDL model as per Section 4.3
spyall.new <- dLagM::ardlDlm(formula = SPY ~ TYS + VIX + IJC, 
                             data = data.frame(merge(returns_spy, diff_tys, weekly_vix[-1,], ijc_growth)["/2022-05",]), 
                             p = 4, q = 4, 
                             remove = list(p = list(TYS = c(1:4)))) # Remove all lags of TYS

stargazer(spyall.new$model, type = "text",
          title = "ARDL Regression Using dLagM Package", 
          dep.var.labels.include = F, 
          column.labels = "SPY.t")

## 
## ARDL Regression Using dLagM Package
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                SPY.t           
## -----------------------------------------------
## TYS.t                        0.014***          
##                               (0.005)          
##                                                
## VIX.t                        -0.006***         
##                              (0.0001)          
##                                                
## VIX.1                        0.005***          
##                              (0.0002)          
##                                                
## VIX.2                         0.00003          
##                              (0.0003)          
##                                                
## VIX.3                         -0.0004          
##                              (0.0003)          
##                                                
## VIX.4                        0.001***          
##                              (0.0002)          
##                                                
## IJC.t                        0.010***          
##                               (0.001)          
##                                                
## IJC.1                        -0.011***         
##                               (0.001)          
##                                                
## IJC.2                        0.007***          
##                               (0.001)          
##                                                
## IJC.3                          0.001           
##                               (0.001)          
##                                                
## IJC.4                        -0.004***         
##                               (0.001)          
##                                                
## SPY.1                         -0.021           
##                               (0.030)          
##                                                
## SPY.2                         -0.050*          
##                               (0.030)          
##                                                
## SPY.3                        -0.111***         
##                               (0.029)          
##                                                
## SPY.4                          0.026           
##                               (0.018)          
##                                                
## Constant                     0.008***          
##                               (0.001)          
##                                                
## -----------------------------------------------
## Observations                   1,165           
## R2                             0.679           
## Adjusted R2                    0.674           
## Residual Std. Error      0.014 (df = 1149)     
## F Statistic         161.666*** (df = 15; 1149) 
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

5.2.2 Forecast ARDL Model

# Create matrix containing values of exogenous variables during forecast period
# Number of columns = forecast period, Number of rows = exogenous variables
x.new = rbind(as.vector(diff_tys["2022-06",]),
              as.vector(weekly_vix["2022-06",]),
              as.vector(ijc_growth["2022-06",]))

f_ardl <- dLagM::forecast(model = spyall.new, x = x.new, h = 4, interval = T, level = 0.95)

stargazer(f_ardl$forecasts, 
          type = "text", 
          title = "Forecasted SPY Returns for June 2022 Using ARDL Model",
          summary = F)

## 
## Forecasted SPY Returns for June 2022 Using ARDL Model
## ========================
##   95% LB Forecast 95% UB
## ------------------------
## 1 -0.047  -0.018  0.009 
## 2 -0.056  -0.029  -0.002
## 3 -0.004  0.025   0.052 
## 4 -0.037  -0.010  0.019 
## ------------------------

To calculate the the forecast accuracy measures using the accuracy function, we need to “trick” the function into thinking it is an object of class “forecast”.

f_ardl.new <- structure(list(level = 95,
                             mean = f_ardl$forecasts$Forecast,
                             lower = f_ardl$forecasts$`95% LB`,
                             upper = f_ardl$forecasts$`95% UB`,
                             x = returns_spy["/2022-05",],
                             fitted = c(rep(NA, 4), spyall.new$model$fitted.values),
                             residuals = spyall.new$model$residuals),
                        class = "forecast")

stargazer(forecast::accuracy(object = f_ardl.new, x = returns_spy["2022-06"]), 
          type = "text", 
          title = "Evaluation of ARDL Model Forecast")

## 
## Evaluation of ARDL Model Forecast
## ===================================================
##                ME   RMSE   MAE   MPE    MAPE  MASE 
## ---------------------------------------------------
## Training set   -0   0.014 0.010         Inf   0.578
## Test set     -0.012 0.033 0.032 62.250 62.250 1.878
## ---------------------------------------------------

Plot forecasted and actual SPY at level (price data):

f_ardl.level <- zoo(x = cbind(fspy = as.vector(last(spy_close["2022-05",])) * cumprod(1 + f_ardl.new$mean),
                              fspy_lower = as.vector(last(spy_close["2022-05",])) * cumprod(1 + f_ardl.new$lower),
                              fspy_upper = as.vector(last(spy_close["2022-05",])) * cumprod(1 + f_ardl.new$upper)),
                    order.by = index(spy_close["2022-06",]))

plot.zoo(cbind(f_ardl.level, spy_close["2022-06",], spy_close["2022-01/2022-05",]), 
         plot.type = "single", 
         col = c("red", "gray", "gray", "blue", "black"), lwd = 2, 
         ylab = "SPY Weekly Close Price", main = "Actual and Forecasted SPY Price Using ARDL Model")

legend(x = "bottomleft", 
       legend = c("Forecasted Price", "Lower and Upper Bound", "Actual Price June 2022", "Actual Price in Train Period"), 
       col = c("red", "gray", "blue", "black"), lwd = 1.5)

The forecast assumes that we have data of the exogenous variables in the forecast period. In reality, we do not have such data, and we would require forecasts of the exogenous variables, which would also decrease the accuracy of the forecasted points and the prediction intervals. This highlights the difficulty of producing accurate forecasts of short-term stock market movements, especially for longer periods.

6 Final Remarks

The ARDL model created in Section 4.3 suggested that past values of the 10Y/2Y Treasury Yield Spread does not explain current SPY returns. This is surprising as I had thought that that knowing the past few yield spread data would help in predicting the current SPY returns.

The models in this project does not consider the effects of structural breaks, where relationships between variables may change due to exogenous factors. A regression model using rolling-window estimation is one of the simpler methods to capture these changes. I would explore the effectiveness of rolling-window estimation in another project and see how it performs compared to a fixed-window estimation.

References

Chicago Board Options Exchange, CBOE Volatility Index: VIX [VIXCLS], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/VIXCLS, July 16, 2022.

Federal Reserve Bank of St. Louis, 10-Year Treasury Constant Maturity Minus 2-Year Treasury Constant Maturity [T10Y2Y], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/T10Y2Y, July 15, 2022.

Kenton, W. (2022). The S&P 500 Index: Standard & Poor’s 500 Index. Retrieved from https://www.investopedia.com/terms/s/sp500.asp

U.S. Employment and Training Administration, Initial Claims [ICSA], retrieved from FRED, Federal Reserve Bank of St. Louis; https://fred.stlouisfed.org/series/ICSA, July 15, 2022.

Walden, S. (2022). Inverted Yield Curve. Retrieved from https://www.investopedia.com/terms/i/invertedyieldcurve.asp.

Effects of Bond Yields, VIX and Initial Jobless Claims on S&P 500

Tan Zheng Liang

2022-07-21

1 Introduction

2 Packages Required

3 Description of Variables

3.1 SPDR S&P 500 ETF

3.2 10Y/2Y Treasury Yield Spread

3.3 CBOE Volatility Index

3.4 U.S. Initial Jobless Claims

3.5 Summary

4 Building the Regression Models

4.1 Autoregressive Moving Average (ARMA) Model

4.2 Autoregressive Distributed Lag (ARDL) Model With One Independent Variable

4.2.1 Regress SPY Returns on Treasury Yield Spread

4.2.2 Regress SPY Returns on VIX

4.2.3 Regress SPY Returns on Initial Jobless Claims

4.3 ARDL Model With All Independent Variables

5 Model Evaluation

5.1 ARMA Model

5.2 ARDL Model

5.2.1 Re-estimating ARDL Model

5.2.2 Forecast ARDL Model

6 Final Remarks

References